
#####################################################
#
#               Code to replicate 
#  Affluence and influence in a social democracy 
#   
#    Ruben Berge Mathisen (University of Bergen)
#
#       Tables and Figures for Manuscript
#                February 2022
#
#   Original analyses conducted in R version 3.5.1 
#
#####################################################



#####################################################
######## Packages and data prep #####################
#####################################################

library(car)
library(dplyr)
library(tidyr)
library(readxl)
library(kableExtra)
library(ggplot2)
library(ggrepel)
library(gridExtra)
library(stringr)
library(stargazer)
library(cowplot)
library(arsenal)
d <- read_excel("Data for Affluence and influence in a social democracy Feb 2022.xlsx")

# Logit-transform variables
logitTransform <- function(p) { log(p/(1-p)) }
d$all_logit <- logitTransform(d$all)
d$inc_p90_logit <- logitTransform(d$inc_p90)
d$inc_p70_logit <- logitTransform(d$inc_p70)
d$inc_p50_logit <- logitTransform(d$inc_p50)
d$inc_p30_logit <- logitTransform(d$inc_p30)
d$inc_p10_logit <- logitTransform(d$inc_p10)
rec <- function(x) {
  x <- ifelse(x>=1,0.9999,ifelse(x<=0,0.0001,x))
  logitTransform(x)  
}
d <- d %>% mutate_at(vars(pred_i90e90:pred_i10e10),funs(rec))

# Remove half-adopted policies 
d <- d %>% filter(outcome!=0.5)

# Remove constitutional issues 
d <- d %>% filter(constitutional_change!="yes")

# Remove items from academic surveys 
d <- d %>% filter(academic_survey==0)


#####################################################
######## Functions ##################################
#####################################################

# Functions for calulcating table and figure quantities
ests <- NULL
eff<-function(x){
  m<-glm(outcome~x,data=d1,family="binomial")
  ests <- data.frame(rbind(ests,round(summary(m)$coefficients[2,1],2)))
  p <- summary(m)$coefficients[2,4]
  se <- round(summary(m)$coefficients[2,2],2)
  p1 <- car::recode(p,"lo:0.009999='***'; 0.01:0.049999='**'; 0.05:0.099999='*'; 0.1:hi=''")
  result <- paste0(ests, p1)
} 
eff_se<-function(x){
  m<-glm(outcome~x,data=d1,family="binomial")
  ests <- data.frame(rbind(ests,round(summary(m)$coefficients[2,1],2)))
  se <- round(summary(m)$coefficients[2,2],2)
  result <- paste0(ests," (",se,")")
} 
eff_num<-function(x){
  m<-glm(outcome~x,data=d1,family="binomial")
  ests <- data.frame(rbind(ests,round(summary(m)$coefficients[2,1],2)))
  ests
} 
se<-function(x){
  m<-glm(outcome~x,data=d1,family="binomial")
  se <- round(summary(m)$coefficients[2,2],2)
  se
} 
p<-function(x){
  m<-glm(outcome~x,data=d1,family="binomial")
  p <- summary(m)$coefficients[2,4]
  substring(as.character(ifelse(p<0.001,"<<.001", round(p,3))),2)
} 
sup20<-function(x){
  m<-glm(outcome~x,data=d1,family="binomial")
  round(predict(m,newdata=data.frame(x=logitTransform(0.2)),type="response"),2)
} 
sup80<-function(x){
  m<-glm(outcome~x,data=d1,family="binomial")
  round(predict(m,newdata=data.frame(x=logitTransform(0.8)),type="response"),2)
} 
fact<-function(x){
  m<-glm(outcome~x,data=d1,family="binomial")
  diff<- predict(m,newdata=data.frame(x=logitTransform(0.8)),type="response")/predict(m,newdata=data.frame(x=logitTransform(0.2)),type="response")
  round(diff,1)
} 
n<-function(x){
  m<-glm(outcome~x,data=d1,family="binomial")
  length(m$residuals)
} 


#####################################################
##### Tables and Figures (same order as in paper) ###
#####################################################

# Table 1
d1 <- d
tab<-data.frame(all=t(d1 %>% summarise_at(vars(all_logit),funs(eff,se,sup20,sup80,fact,n))),
                inc10=t(d1 %>% summarise_at(vars(inc_p10_logit),funs(eff,se,sup20,sup80,fact,n))),
                inc30=t(d1 %>% summarise_at(vars(inc_p30_logit),funs(eff,se,sup20,sup80,fact,n))),
                inc50=t(d1 %>% summarise_at(vars(inc_p50_logit),funs(eff,se,sup20,sup80,fact,n))),
                inc70=t(d1 %>% summarise_at(vars(inc_p70_logit),funs(eff,se,sup20,sup80,fact,n))),
                inc90=t(d1 %>% summarise_at(vars(inc_p90_logit),funs(eff,se,sup20,sup80,fact,n))))
tab1 <- data.frame(t(tab))
row.names(tab1) <- c("All","P10","P30","P50","P70","P90")
kable(tab1, align="c",booktabs=T,linesep = '\\addlinespace', 
      col.names = c("Effect (logit coefficient)",
                    "Standard error",
                    "Predicted probability of policy change at 20% support",
                    "Predicted probability of policy change at 80% support",
                    "Relative change in probability",
                    "N")) %>%
  column_spec(1,width="0.9in") %>%
  column_spec(2:3,width="0.6in") %>% 
  column_spec(4:5,width="0.8in") %>% 
  column_spec(6,width="0.5in") %>% 
  column_spec(7,width="0.4in") %>% 
  pack_rows("Income percentile",2,6) %>%
  save_kable("A&ISD Table 1.doc")


# Table 2
d1 <- d %>% filter(diff_inc90_inc50>0.1)
tab1<-data.frame(inc50=t(d1 %>% summarise_at(vars(inc_p50_logit),funs(eff,se,sup20,sup80,fact,n))),
                 inc90a=t(d1 %>% summarise_at(vars(inc_p90_logit),funs(eff,se,sup20,sup80,fact,n))))
d1 <- d %>% filter(diff_inc90_inc10>0.1)
tab2<-data.frame(inc10=t(d1 %>% summarise_at(vars(inc_p10_logit),funs(eff,se,sup20,sup80,fact,n))),
                 inc90b=t(d1 %>% summarise_at(vars(inc_p90_logit),funs(eff,se,sup20,sup80,fact,n))))
tab <- cbind(tab2,tab1)
tab1 <- data.frame(t(tab))
row.names(tab1) <- c("Poor (P10)",
                     "Affluent (P90) ",
                     "Middle (P50)",
                     "Affluent (P90)")
kable(tab1,type="word", booktabs=T, align="c",linesep = '\\addlinespace', col.names = c("Effect (logit coefficient)","Standard error", "Predicted probability of policy change at 20% support","Predicted probability of policy change at 80% support","Relative change in probability","N")) %>%
  column_spec(1,width="1.4in") %>%
  column_spec(2,width="0.5in") %>% 
  column_spec(3,width="0.55in") %>% 
  column_spec(4:5,width="0.8in") %>% 
  column_spec(6,width="0.5in") %>% 
  column_spec(7,width="0.4in") %>% 
  pack_rows("Poor vs. affluent",1,2) %>%
  pack_rows("Middle vs. affluent",3,4) %>%
  kable_styling()  %>%
  save_kable("A&ISD Table 2.doc")

# Figure 1
diff <- d %>% 
  filter(diff_inc90_inc50>0.1)
m<-glm(outcome~inc_p90_logit, data=diff,family = "binomial")
newdata <- with(diff, data.frame(inc_p90_logit= logitTransform(seq(0.05,0.95,by=0.01))))
tab1 <- data.frame(pred=predict(m,newdata,type = "response", se.fit=T), labels= seq(0.05,0.95,by=0.01))
tab1$group <- "Affluent"
m<-glm(outcome~inc_p50_logit, data=diff,family = "binomial")
newdata <- with(diff, data.frame(inc_p50_logit= logitTransform(seq(0.05,0.95,by=0.01))))
tab2 <- data.frame(pred=predict(m,newdata,type = "response", se.fit=T), labels= seq(0.05,0.95,by=0.01))
tab2$group <- "Middle"
tab3 <- rbind(tab1,tab2)
tab3$group <- factor(tab3$group)
tab3$group <- factor(tab3$group,levels(tab3$group)[c(1,2)])
fig1<-ggplot(tab3,aes(x=labels,y=pred.fit,linetype=group)) +
  geom_line(size=0.8) +
  # geom_ribbon(aes(ymin=pred.lwr,ymax=pred.upr),alpha=0.07) +
  scale_y_continuous(breaks=seq(0,1,0.2), labels=function(mean) paste0(mean*100,"%")) +
  scale_x_continuous(limits=c(0,1), breaks=seq(0,1,0.2), labels=function(labels) paste0(labels*100,"%")) +
  labs(x="Percent favoring policy change",y="Predicted probability of policy change", 
       linetype=NULL, 
       color=NULL,
       title=NULL) +
  theme_bw() +
  scale_linetype_manual(values=c("solid","longdash")) +
  coord_cartesian(ylim=c(0,1)) +
  theme(legend.position = "None",
        legend.background = element_rect(size=0.2, colour ="black"),
        text = element_text(size = 14),
        panel.grid = element_blank()) +
  guides(linetype = guide_legend(override.aes = list(size = 0.4))) +
  annotate("text",0.5,0.6,label="High income (90th percentile)",size=3.5) +
  annotate("text",0.7,0.15,label="Median income (50th percentile)",size=3.5)
ggsave(file="A&ISD Figure 1.pdf",
       plot = fig1,
       width = 4.8, height = 4.3)

# Figure 2
d1 <- d %>% filter(diff_inc90_inc10>0.08&topic=="Economy")
tab1<- data.frame(lab=t(d1 %>% summarise_at(vars(inc_p10_logit,inc_p90_logit),funs(eff_se))),
                  eff=t(d1 %>% summarise_at(vars(inc_p10_logit,inc_p90_logit),funs(eff_num))),
                  topic=paste0("Economic policy\n(n=",nrow(d1),")"))
d1 <- d %>% filter(diff_inc90_inc10>0.08&topic=="Moral")
tab2<- data.frame(lab=t(d1 %>% summarise_at(vars(inc_p10_logit,inc_p90_logit),funs(eff_se))),
                  eff=t(d1 %>% summarise_at(vars(inc_p10_logit,inc_p90_logit),funs(eff_num))),
                  topic=paste0("Moral policy\n(n=",nrow(d1),")"))
d1 <- d %>% filter(diff_inc90_inc10>0.04&topic=="Foreign/Security")
tab3<- data.frame(lab=t(d1 %>% summarise_at(vars(inc_p10_logit,inc_p90_logit),funs(eff_se))),
                  eff=t(d1 %>% summarise_at(vars(inc_p10_logit,inc_p90_logit),funs(eff_num))),
                  topic=paste0("Foreign/Security policy\n(n=",nrow(d1),")"))
d1 <- d %>% filter(diff_inc90_inc10>0.04&topic=="Other")
tab4<- data.frame(lab=t(d1 %>% summarise_at(vars(inc_p10_logit,inc_p90_logit),funs(eff_se))),
                  eff=t(d1 %>% summarise_at(vars(inc_p10_logit,inc_p90_logit),funs(eff_num))),
                  topic=paste0("Other\n(n=",nrow(d1),")"))
tab <- rbind(tab1,tab2,tab3,tab4)
tab$percentile <- ifelse(grepl("inc_p10",rownames(tab)),"Low income (10th percentile)","High income (90th percentile)")
tab$percentile <- factor(tab$percentile,levels(factor(tab$percentile))[c(2,1)])
fig2<-ggplot(tab,aes(topic,eff,fill=percentile)) +
  geom_bar(width=0.75,stat="identity",position="dodge",color="black") +
  scale_fill_manual(values=c("grey90","grey30")) +
  geom_text(aes(label=lab),position=position_dodge(width=0.8),
            size=2.7,
            vjust=ifelse(tab$eff==-0.25,-5.8,-1.3)) +
  geom_text(aes(label=ifelse(percentile=="High income (90th percentile)",as.character(topic),NA),y=1.4)) +
  scale_y_continuous(limits=c(-0.3,1.7)) +
  theme_bw() + 
  theme(panel.grid = element_blank(),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        legend.position = "bottom") +
  labs(x=NULL,y="Opinion-policy link (logistic coefficient)",fill=NULL) 
ggsave(file="A&ISD Figure 2.pdf",
       plot = fig2,
       width = 7.2, height = 4)

# Figure 3
d1 <- d %>% filter((diff_inc90_inc10>0.1 | diff_edu90_edu10>0.1))
norway<-data.frame(coef=t(d1 %>% summarise_at(vars(pred_i90e10,pred_i10e90),funs(eff_num))),
                   se=t(d1 %>% summarise_at(vars(pred_i90e10,pred_i10e90),funs(se))))
# US estimates from Table A3.4 in Gilens (2012, 259): 
# Income P90 and Education P10: 0.41 (0.08)
# Income P10 and Education P90: 0.27 (0.08) 
d1 <- data.frame(country=c(rep("United States",2),rep("Norway",2)),
                 edu=c("High income, low education","High education, low income","High income, low education","High education, low income"),
                 coef=c(0.41,0.27,c(norway$coef)),
                 se=c(0.08,0.08,c(norway$se)))
d1$lab=paste0(d1$coef," (",d1$se,")")
d1$edu <- factor(d1$edu,levels(factor(d1$edu))[c(2,1)])
fig3<-ggplot(d1,aes(country,coef,fill=edu)) +
  geom_bar(stat="identity",position="dodge",color="black",width = 0.7) +
  scale_fill_manual(values=c("grey30","grey90")) +
  geom_text(aes(label=lab),position=position_dodge(width=0.75),vjust=-1.3,size=3.2) +
  labs(x=NULL,y="Opinion-policy link (logistic coefficient)",fill="Subgroup") +
  theme_bw() +
  theme(panel.grid = element_blank(),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank()) +
  scale_y_continuous(limits=c(0,0.7)) +
  geom_text(aes(label=ifelse(edu=="High income, low education",as.character(country),NA),y=0.6))
ggsave(file="A&ISD Figure 3.pdf",
       plot = fig3,
       width = 6.8, height = 3)

###################### END #############################


